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INTRODUCTION 


This  final  report  on  global  weather  modeling  with  vector  spherical 
harmonics  supported  by  the  Air  Force  Office  of  Scientific  Research 
Contract  No.^44620-74-C-0036  under  ARPA  order  2633  describes  a 
third  dimensional  spectral  model  based  on  the  use  of  vector  spherical 
harmonics  and  scalar  spherical  harmonics  in  the  horizontal  and  Legendre 
Polynomials  in  the  vertical  direction.  Thus  at  each  time  step  no  spatial 
derivatives  but  only  time-derivatives  need  to  be  computed.  The  vector 
spherical  harmonics  have  two  distinct  advantages:  no  horizontal  deriva- 
tives appear  at  any  point  of  the  formulation  and  the  functions  have  a 
clearly  defined  direction  at  the  poles.  The  model  implemented  is  what 
corresponds  to  and  eight  layer  version  of  Gates' Mintz-Arakawa  model. 

The  first  chapter  is  a review  of  the  mathematics  which  has  been  an 
attempt  at  a formulation  which  is  new  in  the  sence  that  it  tries  to 
make  the  presentation  as  simple  as  possible  by  not  starting  with  the 
usual  group  theoretical  approach.  The  chapter  on  the  program  design 
has  been  included  for  completeness,  although  it  was  not  supported  by 
this  contract  but  thru  several  IRAD  (Independent  Research  and  Development) 
tasks  by  IBM.  This  includes  both  design  and  implementation  of  the 
general  program  system.  The  programming  done  under  the  contract  has 
been  for  specific  interface  programs  to  enable  us  to  make  computer 
runs  with  real  data. 
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MATHEMATICAL  BACKGROUND 


Since  the  mathematical  formulation  used  in  this  paper  is  not  commonly  used 
in  meteorology,  it  is  useful  to  include  a discussion  of  the  mathematical  tools 
used. 

There  are  several  reasons  why  it  is  natural  to  express  scalar  fields  such 
as  pressure,  temperature  and  so  on  in  spherical  coordinates  in  terms  of 
spherical  harmonics.  The  most  compelling  reasons  are  group  theoretical 
and  we  will  not  discuss  them  here.  The  spherical  harmonics  Y™  ( 6, cp) 
form  a complete  orthonormal  set  of  functions  in  the  co-latitude  6 and  longi- 
tude cp  . This  means  that  any  scalar  function,  i.  e.  , one  where  one  number 
is  associated  with  the  physical  variable  at  each  point  in  space,  can  be  ex- 
pressed as  a series  in  terms  of  the  spherical  harmonics,  Consider,  for 
example,  the  pressure  p to  be  a function  of  6 and  cp  , Then  one  can  write 

p(0,cp)  = Tj  pnl  Y™  (0 , cp)  (1)  • 

m,  i 

where  p™  are  constant  coefficients,  not  to  be  confused  with  the  Legendre 
functions  introduced  later. 

The  functions  Y1]?  (B,cp)  are  defined  in  the  following  way: 

Y^e.cp)  = (-i)m 


2 & + 1 (n  - m) ! 
4 tt  {t  + m) ! 


>£(cos6)eimC*  (2) 


It  should  be  pointed  out  that  the  spherical  harmonics  and  the  Legendre 
functions  are  defined  in  slightly  different  ways  in  different  papers  anu 
that  it  is  important  to  adhere  to  the  definitions  used  here  for  the  subse- 
quent equations  to  be  correct. 


it  takes  on  integral  positive  values  or  zero. 


-L  = 0,  1 , 2 , • • • . For  a given  value  of  b , m can  have  2 1+1  values 

0 + 1 +2  • • • + ' 
u 1 1 _ c i ^ * 


P™  (cos  0)  is  an  associated  Legendre  function  and  has  the  form 


P1)1  (cos  6) 


■ m n d 
sin  0 


m 


d(cos  0) 


[P  (cos  6)] 


m £, 


(3) 


Where  P^  (cos  0)  is  a Legendre  Polynomial  defined  by 


P (cos  9)  = 


d ’ [c os^  0 - 1 1 


2^  .L ! d ( cos  0) 1 


(4) 


The  first  few  Legendre  Polynomials  are 


PQ  (cos  0)  = 1 


Pj  (cos  0)  = cos  0 


2 

n / 3 COS  0-1 

P^  (cos  0 ) = ^ 


P^  (cos  0 ) = 


5 cos  G-3  cos  0 


2 


r 


^ WWWHu  M IIJ A!-M  IWW  IP  ■ 


The  first  few  associated  Legendre  functions  arc: 


PJ  <cos  e)  = sin  9 , P*  (cos  0)  = 3 cos  6 sin  0 


PZ  (cos  6)  = 3 sin  0 , pj  = ~ sin  0 (5  cos  0 - ] ) 


P3  (cos  0 ) = 15  sin  0 cos  0 , P3  = 1 5 sin3 


P[  = P. 

I .L 


The  spherical  harmonics  Y ^(0,cp)  satisfy  the  following  orthogonality 
conditions. 


Z TT  TT  1 

rm* 


p p m 5-c  rn 

J J Y ,)  (0.cp)  Y (0  ,cp)  sin  0 d 0 d cp  = 6 ,6  , 

oo  ^ l 


The  symbol  6^  ml  is  a Kronicker  delta  and  is  equal  to  zero  if  m / m* 
and  equal  to  one  if  m = m1  . 
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lhe  functions  (0(Cp)  can  then  be  lormed  from  the  examples  above  by 

the  use  of  equation  (2).  The  spherical  harmonics  satisfy  the  equation 

v2  Y-fe.cp,  =-— UJ_!)  v?(e.v) 

r 

U a scalar  field  is  a function  of  the  three  spherical  coordinates  r , 0 , cp  , 
and  of  time  it  can  be  expanded  in  the  following  form 

P( r , 0 , cp , t)  = £ p™  (r,  t)  Ym  (6  ,cp)  (6) 

m,i  ' x ; 

Now  the  expansion  coefficients  are  functions  of  r and  t 


The  gradient  of  the  pressure  is  a vector  field  so  it  is  appropriate  to  look 
for  a representation  of  a vector  field  in  terms  of  functions  of  0 and  cp 

which  are  the  logical  counterparts  of  the  scalar  spherical  harmonics 

yT  (e,cp)  . 

As 

Let  us  form  the  gradient  of  a typical  term  in  the  expansion  for  the  pres- 
sure. We  will  for  the  time  being  omit  the  time  dependence  which  has  no 

influence  on  the  contents  of  the  next  few  pages.  Consider  a term  of  the 
form 


i j-  .iijiiia.  iii  y J.  a.|A".J.‘.i.fl'i"w,i,i!Mi  iy.h  ,,  , ,1.1 , ,. ir  i W WJJP  'tvyrn*!*' 


Ihe  functions  (<^, cp)  can  then  be  formed  from  the  examples  above  by 

the  use  of  equation  (2).  The  spherical  harmonics  satisfy  the  equation 


rm 

jl 


( e- 9) 


1) 


Y?1 


(O.cp) 


If  a scalar  field  is  a function  of  the  three  spherical  coordinates  r , 0 , cp  , 
and  of  time  it  can  be  expanded  in  the  following  form 


P(r,  e.cp,  t)  = L P™  (r,  t)  Y,f"  (G  , cp) 
m ,1  L 


Now  the  expansion  coefficients  are  functions  of  r and  t . 


The  gradient  of  the  pressure  is  a vector  field  so  it  is  appropriate  to  look 
for  a representation  of  a vector  field  in  terms  of  functions  of  6 and  cp 
which  are  the  logical  counterparts  of  the  scalar  spherical  harmonics 
Y™  (0.cp)  . 

.A/ 

Let  us  form  the  gradient  of  a typical  term  in  the  expansion  for  the  pres- 
sure. We  will  for  the  Lime  being  omit  the  time  dependence  which  has  no 
influence  on  the  contents  of  the  next  few  pages.  Consider  a term  of  the 
form 
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.JP-i  M.  UUU-W. 


(5) 


(6) 


■ ,.,.1 


The  gradient  in  spherical  coordinates  is  formed  by  applyi 


ng  tho  operator 


A A 

0 0 


v = £ 0 + _«  _iL  + _j£L 

i'  dr  r 3 B 


r sin  0 3 cp 


(7) 


A A 


A 


Where  er  , efe  , and  e^  are  an  orthogonal  set  of  unit  vectors  pointing 
in  the  direction  of  increasing  r , & , and  cp  . 


First,  consider  only  the  term  Y™  (8,cp) 


V Y™  (b.cp)  = i 


dY 


m 

i + 


3Y™ 


L 9 ^ ® cp  s in  0 3 cp 


(8) 


The  result  is  a vector  field  that  has  0 and  cp  components,  i.e.  , which  is 
tangential  to  the  surface  of  a sphere.  We  can  then  use  this  as  a part  of  a 
general  vector  field.  Let  us  define  a field  B^  (0,cp)by 


(0  .cp)  = 


V Y 


M 


1 


V"~L(L  + 1) 

The  reason  for  the  factor 


3YM 

A Y L ^ A 
c„  ts-  + e 


1 


3 Y 


M 


L VL(L  + 1)  L~e  d 9 > sin  0 9cP 

1 


(9) 


will  be  discussed  later. 


7L(L  + 1) 


We  will  use  the  notation 


YL  = %/L(L  + l)  (10) 

M 

® »V)  WH1  now  be  one  of  our  basic  fields.  One  can  see  from  equations 

(2),  (3)  and  (9)  that  if  M = 0 , (0  ,cp)  is  only  a function  0 and  only  has  a 

component  in  the  0 direction. 

L(0)  thus  has  the  same  value  on  any  point 
on  a longitude  circle  and  everywhere  points  in  a direction  parallel  to  the 
lines  of  constant  longitude. 


M 


Now  consider  the  gradient  of  f(r)  Y /(S.cp) 


7 f (r)  Y^(0(Cp)  = er  Y^(e,cp)  + f(r)  7 Y^(e>cp) 


d f ( r ) A VM 

dr 


er  Yi><»  + r 


L 


M 


f(r)  B"  (e.cp) 


(11) 


We  will  now  define  another  vector  field  A^(p,cp)  = e Y^(&,cp) 


(12) 


M 

A L (0,cp)  is  a vector  field  which  points  in  the  radial  direction.  We  can  then 
rewrite  (12)  in  the  form 

Y 


?f(r)  y“(0.cp)  - A“(e.p)  + -i-  f (r)  b“(0  ,<p) 

A L (0,cp)  and  B^(0,cp)  are  not  sufficient  to  form  the  basis  for  the 
representation  of  a general  vector  field  which  is  a function  of  0 and  cp  , 
simply  because  vector  fields  are  not  all  expressible  as  gradients  of 


(13) 


<0) 


scalars.  This  can  also  be  seen  by  looking  at  the  special  case  A®  and 

By  • -As  we  have  seen,  A T points  in  the  radial  direction  and  B° 

L points 

in  the  north  south  direction,  so  in  general,  we  in  addition  need  an  east- 
west  component. 


This  leads  us  to  look  for  an  additional  field  which  is  perpendicular  to 

both  and  . The  vector  product  of  two  vectors  is  perpendicular 
, . M . . -A 

to  each.  A ^ is  proportional  to  e^_  so  let  us  define  a new  vector  field 


C^(9,cp)  by 


(6,cp)  - 


•i  er  X B^  (0,cp) 


( 14) 


(The  reason  for  the  factor  i will  become  apparent  later.  ) Now, 


A 

e 


c: (0,cp)  = -i~  x 

IT 


_r 


1 


9Y 


A A _ 

60  9e  ecp  sin  0 9 0 


M 


„ A A 
r * % = e 


A A A 

e X e = -eQ 
r cp  0 


Thus 


c1^  (e, cp) 


9 Y 


M 


cp  90 


1 


9 Y 


M 


0 sin9  ^cp 


( 15) 


( 16) 


( 17) 


We  see  that  only  has  a cp  component.  We  will  show  later  that 
the  set  of  A^(6,cp)  , B^(6,cp)  and  C^(0,cp)  form  a complete 
or thonoi  mal  set  of  vector  functions  in  8 and  cp  and  consequently  any 
mathematically  well-behaved  vector  field  which  is  a function  of  the  indepen- 
dent variables  0 and  cp  can  be  expressed  as  a series  of  these  functions.  As 
an  example,  take  the  velocity  field  V of  a fluid.  Let  V be  a function  of  the 
spherical  coordinates  r , 9 and  cp  and  of  time  t . Then  V(r,0,cp,t)  can 
be  expressed  as 


V(r , 9,  cp  , t) 


\ M . . M /Q  . 

a (r,  t)  A . ( 0 , cp ) 

M,  L 


(r , t) 


(G,cp  ) 


+ Z.  C l!  (r*  C L (6,Cp) 

M,  L 


Now  that  we  have  found  a representation  that  can  be  related  in  a simple 
and  logical  way  to  the  gradient  of  a scalar  field,  we  must  test  its  utility. 

In  order  to  do  this  we  will  examine  the  curl  and  divergence  operations  in 
terms  of  our  representation.  Since  the  curl  of  a vector  field  is  itself  a 
vector  field,  the  result  of  the  curl  operation  should  be  in  terms  of  the 
•A  l > B and  C ^ functions.  Since  the  divergence  of  a vector  field  is 
a scalar  field,  the  divergence  of  A1^  , B^  and  should  be  in  terms 

of  scalar  spherical  harmonics  Y ^ . It  is  not  clear  at  this  point  how 
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simple  the  results  will  be  in  comparison  with  the  conventional  way  of 
doing  these  operations  in  spherical  coordinates  where  all  the  deriva- 
tives with  respect  to  the  independent  variables  occur.  In  addition,  the 
conventional  formulation  contains  numerous  trigonometric  expressions. 


We  begin  with  the  divergence. 


Consider 


f(r)  (6,cp)  = f(r)  e^  (&,cp) 


V • f(r)  A l (9.CP)  = 


; _d_ 

'r  dr 


r sin  6 dcp 


' e f(r)  Y^(e.cp)  (19) 


One  must  note  that  there  is  an  important  difference  between  Cartesian  and 
spherical  coordinates.  In  Cartesian  coordinates  the  derivatives  of  the 

AAA 

unit  vectors  1 , j , k are  zero.  In  spherical  coordinates  the  derivatives 

of  e , e 0 and  e are  in  general  not  zero  since  e , en  and  e change 

r 0 cp  ° 

direction  from  point  to  point.  It  .is  customary  in  meteorological  work 
when  going  to  spherical  coordinates  to  write 


A ‘ = 1 9 

d x r sin  0 ^ cp 


This  is  a dangerous  practice  since  it  ignores  the  above-mentioned  fact. 


We  will  here  list  the  derivatives  for  future  reference: 


p 


I I.  Mil  .!•  II  Mil..*--.  ■ < ^ 

r 

i 

F 


li  we  break  the  operator  V into  two  parts 


AS,! 
e — + — 

r S r r 


A 

A _S_  ( cp  S 

C0  S0  sin  G Sep 


V1  + V2 


whe  re 


A 

e — — and  V. 
r o r , 


i A S 

r °0  se 


sin  0 Bg 


We  find  V.  f(r)  J3J,  (9,g) 


<V(r»  * Bl  (Q.<?>  +f(r)v2-  B^(a.cp) 


We  kn  that 


2 M 
V Y (P.  cp) 


\ YT  (S-V)  = 7 / y“(6  ,cp> 


Thus  V2  . Bl 


V_ 

yl  2 


M 

r v Y 
2 L 


iz 


Ik  ym 

r L 


and  consequently , since  VJ  f(r)  = ef  ~~  and  B1^  are  perpendicular 
we  (j;et 


V • £(r)  ( 0, g)  = 


f(r)  Y™  (0fCp) 


(24) 


Lastly,  we  want  to  find  V ■ f(r)  C,  ( 9 , g) 


12 


- . Ji,.. 


X,  U 1.1. 


Since 


v ■ (a  A)  = a v - A + A • y a 


(25) 


V • (f(r)  C^)  = f(r)V  • C1^  + • Vf(r) 


L 


M A 


.M 


But  Vf(r)  - er  and  is  tangential  so  C1^  • Vf(r)  = 


and  V • (f(r)  ) = f(r)  V-  C1^ 


Now  consider 


V • C 


L = -fL*  ’ (rVn!> 


(26) 


Using  the  identity 


V • (A  X B)  = B • V X A - A • V x B 


(27) 


we  get 


V • C 


M 


V YL 


f 

K 


V X r c 


r e • ^ V X 


V Y 


Lj 


(28) 


whe  re 


VXre  = V X r 
r 


0 and  V x ( V a ) = 0 
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nm«ppfpip«|^IWiimV.  J W«(u  J-  i "Wf  , ,1  m j i, .,, p, * |W. 


thus 


and 


(29) 


V-  f(r)  C ^ (6  , cp)  = 0 
We  collect  the  results  here 


V'  f(r)  = D(2)  f(r)  Y^O.cp) 

(20) 

V-  f(r)  B“(eitp)  = £(r)  y«j9i,| 

v 2 4 ) 

V • f(r)  C^(e,cp)  = 0 

(30) 

when 

U<L>  = < Jr  * r > 

(31) 

We  have  now  found  the  divergence  expressions  and  they  do  indeed  meet 
the  criteria  of  simplicity.  They  also  show  the  logical  interrelationship  ' 
between  on  one  hand  and  A^  and  C ^ on  the  other. 

We  have  found  two  very  important  results.  The  gradient  on  a scalar  field 
expressed  in  terms  of  scalar  spherical  harmonics  Y1^  gives  vector 
spherical  harmonics  A^  and  and  the  divergence  of  the  vector  A^'s 
spherical  harmonics,  and  B^'s  give  scalar  spherical  harmonics. 

Furthermore,  the  simplicity  of  these  expressions  is  striking,  since  only 
derivatives  with  respect  to  r occur  and  the  angular  part  of  the  functions 
just  occur  as  multiplying  factors.  No  derivatives  of  the  angular  part 
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- • - — - * 


or  trigonometric  factors  enter  anywhere.  This  makes  it  very  easy  to 
eliminate  the  angular  part  from  equations. 


If  we  can  now  show  something  similar  for  the  curl,  we  will  have  estab- 
lished both  the  connection  between  the  scalar  and  the  vector  functions 
and  between  the  vector  functions  themselves. 


We  start  with 


v X f (r ) A ^ ( e.cp) 


M 


V X f(r)  A ( O.cp)  = 


V X [ f(r)  e Y1^  ] 
r L 


VX(cp.A)  =VcpXA  + cp  V X A 


A = f(r)  e 


_ vM 
tp  - Yl 


V X (f(r)  a“|  = V y'J1  X £ f(r)  + V X (e  f(r)) 

i_  r L.  r 


V X (e  f(r))  - V f(r)  X e + f(r)  V X e 


V f(r)  = e -~— 
r dr 


V X e =0 
r 


Thus 


V x f(r)  A 


M 


Y 


L 


-f(r)er  iVY^  = i f(r ) C**  (6,cp) 


(31) 
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..  tCr  ± 


v X ( A X B ) - A 7.11  - B 7. A*  + ( B • V ) A - { A ■ V ) B 


7 X -i  (o  x B^) 
' r L 


V X (er  X B^)  = cr  V • B^J  - B^  V • er  + (B^  • V ) cf  - (er  • V)  B^ 


A 5 A , M 

e • V = — e • 7)  B J 

r 3r  ' r ' L 


= 0 


V • e 


V • B 


M 


-y  yM 
' L 1 L 


(b^.  ?)£_  -1[^Aym  + a 


a VM 

] a yl 


r r 0 50  L cp  sin  0 5 cp 

M 


j.iLLL  + S _JL._L.jS 

r 0 ?cp  cp  sin  0 Bcp  r 


8 y“ 

[ L 3 


5 Y 


L 5 T A 


5 

1 rA  L , A 1 

Le  ___ — + e 


5 Y 


M 

L 


ftG  3 0 . 2 5 cp  Bcp  r Z ‘■'"0  5 0 

o in  o r 


cp  sin0  5 cp 


L=  B“ 

r L 


Thus 


V x £(r)  C"4  ( 0,  cp)  = iv 


f(r)  aM 


M 


L r 


L ( ©.CP)  + M -JT  +7)  f(r)  B"  ( e .cp)  ) (33) 


We  have  Lhus  derived  the  gradient  divergence  and  curl  expressions.  We  list 
them  here  for  reference. 
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TT  -pr.wTw-r-w^-r^-^.  i i , ' 'u  «•  M'H'vm1'  -.*• 


V f(r)  Y M ( e.cp) 

V • f(r)  A^  (G.cp) 

V • f(r)  B^  (6.CP) 

V • f(r)  (6.CP) 

V X f(r)  A1^  (G.cp) 

V X f(r)  B^  (G.cp) 

V X f(r)  (G.cp) 


•y 

= -f-  f(r)  A^(G.cp)  + £(r)  B^(0.cp) 

' <■£-  + f)  f(r)  Y“  («.<(,) 

v 

■I-*  /■  t \ K4  . p.  v 

= - — i(r)  * L (6,cp) 

= 0 


= - l-~  f(r)  C^(B,g) 

= 1 ( 7Tr~  + r * C L 


f(r)  A^(G.cp) 


+ ~)  f(r)  B1^  (&.CP  ) 


It  is  difficult  to  imagine  a simpler  set  of  expressions  in  spherical  coordi- 
nates than  this.  It  has  been  customary  to  use  a set  of  the  form 


V 


1<  OMi  Y J_  , y i €g 


■V  L 


(i- . 0 


Y l (e.cp) 


A 

e 

cp 


The  reader  is  urged  to  try  the  above  derivations  in  terms  of  these  functions. 
Two  things  happen;  it  is  not  possible  (at  least  in  our  experience)  to  eliminate 
derivatives  in  6 and  various  trigonometric  expressions;  the  functions  do 
not  reproduce  the  way,  for  example,  the  curl  expressions  above  do. 
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The  functions  above;  A™,  b“  and  c“  have  the  advantage  that  it  is  possible  to 
visualize  their  appearance.  They  are,  however,  not  what  one  calls  vector 
spherical  harmonics  at  all  but  linear  combinations  of  these.  So,  while 
they  are  a convenient  way  of  representing  real  data  and  in  a very 
economical  way  express  the  characteristic  features  atmospheric  motion, 
for  example  C±  represent  the  zonal  motion  and  C®'s  (itfO)  represent  Rossby 
waves, they  are  not  easy  to  use  when  one  deals  with  non-linearities.  For 
these  it  is  more  convenient  to  use  the  actual  vector  spherical  harmonics, 
the  derivation  of  which  from  the  a“'s,  B^'s  and  c“'s  we  sketch  below. 


Introducing  the  following  operators 


9 £ - i C-  - 


- ' l C' 

■'  u U 


V 5 


r\  J 

e/  z t 

on 


s _ \ 


-“A 


V*- 


! 

aJ 


’ll 


(34) 


one  can  show  (1)  that 


v -- 1 i yv''  i 

- i l • — • -- r~  *•  - ' r ^ 

L-1  l 1.  t -A  % \(  l.U’V  ' ) J l I 


U.  \ r 


j 


jj.  \ ( \- - V\ )[  H - \)  v'"'  j 


Vl  ' Hi  1-  < ) 


, Lu)N  - 


r, 


IIUh  Mt\  A - ■■  ■ 


■ ? v a v •*»  \ ^ M 


' \ L. 

“■  X 


- -.Vr 


i ( u « ' ^ . \ Sy’ 

'Ml  N 


6 


.6 


C vis*  )\  ■ 

•V  , >X 


~ V'  vu4  r-,  . s*y 

1 , a \ f 1 v.  4f  : _J  U + \ \ rv-  i > 

i 1.  ( i.  V » < ■ *■  / — •* 

u 


\M+>  ly  - r:  \\- 

U+  \ \ c/  * 


v\u h * k ^ \o 


t-  * \£ 

J 


Y 
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if  one  then  introduces  the  following  basis  vectors 


' i a 


t 


1 1 A 


' \ 


^ 0. 


the  del-operator  can  be  written 


V -If  u k,  ■ « 

A * * I 


_ h a 


/■  t'M  ’ £a>  “ ^1'h  */•  *\>l 


m 


Remembering  that  is  defined  as 


£ 


VA 


■V  V 


\u 


'■  Yu. 

and  using  (39)  in  conjunction  with  (35)- (37)  one  finds  after  some 


straight  forward  but  tedious  manipulations  that  can  be  written 


the  the  form 


) ' 


•Hi  , ) 


U- 


■n*  > 


lT,‘ 


Ul. 


Where 


u.  ‘i.  ■ 


) 


A 


l Vi'Nl  "* 1 * 

T « 

>■  - ' - 1 


( u_-  H V l + \i  )i  " y ' ^ 


L^viVw-\)  \ 


v- 1 -i 


and 


- r,  v,  ^ V'  a 
VCLH  < i'\l*  N +\)  V ' - HV  h.\'  l^V.:4.  . | Y p 

V.  0™'.  v-*  * Ml  l -*■  ’’i  \i  ^ ' 1 L )■—  U+' 


(38) 


(39) 


(40) 


(41) 


„m 


LL+1 


\ 

Til-  N yV*  ~ ’ €' 

Ul^  AClUl,  ) J v * ' 


(42) 


20 


WWW* 


If  one  now  repeats  the  above  procedure  for  the  expression  ^ w''}  Y® 
in  - • . (13)  one  finds  that  A™  can  be  bitten 


S’  * £ 


1 . V> 

»-  \ 

i ) u- 


v'-v  - " 


La  i 

^V-.r)V  V- 


(43) 


Remembering  that  i ^')  A™  is  given  by  (31)  one  can  by  using  the  curl 
operator  (39)  on  f(r)  times  (43)  find  C°.  The  result  is 


£ - T?L  = V 

u nu+  \ > J Y > 


fcS L.  . 

lU1^'V'w  'v  c 


. v,  ^ 

v e 


( y-»  frvYA-  ^ *•  ^ 


M -L  A 

U e+1 


(44) 


We  will  call  T™L_^  t°l  an<*  ^LL+1  co^^ect^veiy  T^*  One  can  then  write 


2 V~|J 

t=-\  r 


S-l  « M ^ 


> 


(45) 


when  the  coefficients  h S’  are  defined  by  equations  (41),  (42) 

and  (44)  T™  are  what  is  conventionally  called  vector  spherical  harmonics. 
The  coefficients^  w.  i. ■ u S are  special  cases  oi  what  is  called 
Clebsh-Gordan  coefficients.  Recognizing  this  one  can  then  call  upon 
the  very  powerful  machinery  if  the  three-dimensional  rotation  group  to 
perform  computations.  We  will  here  only  briefly  sketch  how  one  deals 
with  non-linearities . In  addition  to  the  equations  above  one  needs  two 
more.  Equation  (45)  can  be  inverted  and  written 


I 


>W' 


7 h>  \ 


\ u 


(46) 
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-X.:-  ..■■■  ..  ,.■■■^■'■4..^,:.'.  •*»»»** - 


and  the  product  of  two  spherical  harmonics  can  be  written 


Ym  vm’  ..  -Cj1  \ Jill'-*  1 
*1  *1'  “ 1 


V 


' 'll  ■* 


l_ 

: 


h*h 


j i , m 

” '•  ''  c>  v-  «•  v-  ' 1.  *\  • . • V \\+  t;  "x  V 

' \ I 

L 

The  above  equation  and  the  orthonormality  condition  on  the  t“  namely 

41 


(47) 


Vm  «\ 


_ |w|  V y 

* ' i » ® ^ J)u.'  r b J 

^ ^ L 1 frih*  IV  1 ^ 


(48) 


together  with  rules  for  dealing  with  sums  of  products  of  Clebsch-Gordan 
coefficients  form  the  complete  mathematical  machinery  for  applying  vector 
spherical  harmonics  to  meteorology. 


Section  3.  ON  A POSSIBLY  NEW  MODE  OF  ATMOSPHERIC  WAVE  MOTION 


In  order  to  show  the  relative  ease  with  which  the  vector  spherical 
harmonics  can  be  used  to  produce  analytic  results  we  will  do  an  anlaysis 
of  the  conventional  way  of  expressing  two-dimensional  motion  in  terms 
of  a non-divergent  and  an  irrotational  part.  We  might  mention  that  we 
have  tried  to  repeat  this  analysis  using  the  usual  expressions  for 
the  differential  operations  of  vector  calculus  in  spherical  coordinates 
but  not  been  able  to  get  any  results  due  to  the  complexity  of  the 
expressions. 

One  customary  way  of  representing  the  two-dimensional  motion  of  the 
atmosphere  is  to  break  the  velocity  field  up  into  two  parts 

vv  = v X T ^ (1,2) 

A * A 

k is  a unit  vector  in  the  radial  direction  or  in  our  rotation  V - £ r 
^ is  the  two-dimensional  del-operator  V ^ andV  2 are  then  tangential 

-O  -■> 

to  the  surface  of  a spherical  earth.  From  the  definition  s/ ^ andv2 
one  finds 

V * O ^ Vx"  O (3,4) 


We  will  now  look  at  the  problem  in  three  dimensions , and  will  allow 
initially  a vertical  or  radial  component  which  we  eventually  will  suppress. 
Let  us  then  write  down  to  three-dimensional  velocity  fields,  the  first  one 
non-divergent  and  the  second  one  irrotational 


V.  - 1!  cV,0  C(G  ^ (5) 

1 M.u  L K L u L Ku  L L 


with 


V*  V.  - O 


(6) 


and 


u 


v * h <7> 
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f «* 1 ' . '! 1 MW!  M!  yiliflli.il  A w 


with  V,  = & 

1 


(8) 


using  equ.  (20),  (24),  (30),  (31),  (32),  and  (33)  from  Section  1, 
namely 


V-  { '>  = (i > * a ^ ,r)V;  iSlf , 

-iu^OY^ie.y) 

=-0 

^,viA*f6,y  U 

^ --  v ^ i.^M cN^) 

¥>)  = >^L  V'  )AN».y)*;U  -,1 

V Ur  v 


(9) 

(10) 
(11) 


and 


(12) 
(13) 

■)?.%<y)  (w> 

L» 


to  form  (6)  and  (8)  respectively  we  have 


V*7  --  O ^ 


4 


\h 


(15) 


and 


^ \J-'C  -> 


The  radial  part  must  be  zero  all  by  itself  or 


. -X 

\ L.  " 0 


(17) 
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Consequently  the  irrotatlonal  part  contains  no  c“.  In  addition  (16) 
gives: 

^u-  v( ^ + 1 \ (V  ^ 


First  ve  notice  that  since  c“  is  defined  by 

. . ~v 

C v ( b , v ^ V V 9 , w , 

it  corresponds  to  "x'v. 

\4  k 4? 

It  then  follows  that  the  horizontal  part  in  V # represented  by 

X 'A  ) ' t>\  V ) 

is  not  contained  in  the  conventional  way  of  representing  the  divergence- 
less part  of  the  velocity. 

Now  if  we  want  to  impose  the  condition  that  there  is  absolutely  no  vertical 
motion  allowed  equ.  (15)  leads  to  the  conclusion  that 

b“(r,t)  - 0 


Thus  this  condition  eliminates  the  new  field.  In  addition  it  imposes 
on  the  ir rotational  part  i.e. 


i 'V'  's>,Cf  I 


the  condition 


• * V ) > 


J \ V - C 


or  that 


Now,  since  in  reality  there  is  a small  vertical  motion  let  us  allow  this 
and  see  what  conclusions  can  be  made  about  the  magnitudes  of  b“  and 

First  consider  b™,  form  (15)  we  have 


>'  "*J  U v U 
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Let  us  assume  that  a“  is  proportional  to  the  distance  from  the  surface  of 
the  earth  a. 


r H V ' 

Ck  . r ^ 
u c 


V 


I \M 

3^ 


<■1 

c 


TL  AlU 
u 

V 


— L 
V* 


Since  X *\>  \'t  ( the  radius  of  the  earth 

and  we  ignore  the  second  term.  Thus 


v \n 
c>  c, 


A • 


t.  V > ^ V 

V* 


* Vu  c 

a,-  v U <£  V u l 


or 

Now 

and 


V v \* 

'V' 


; * v 

v»>  L 


<\v  \ C C- 


\ w c viUu-«  > ) *0”  u 

Thus  bx  is  one  to  two  orders  of  magnitude  larger  than  a“.  So 

if  the  vertical  motion  represented  by  a“  is  of  the  order  of  cm/s.  bm 

:an  of  order  of  m/s  and  does  not  appear  neglible. 


On  the  other  hand  we  have  from  equation  (18) 


A repetition  of  the  previous  arguments  indicates  that^“  is  an  order 
or  two  of  magnitude  smaller  than^  °. 


Thus  without  having  actual  numbers  fora“  aud  ^ J it  appears  safe  to 
conjecture  that 

v;  >yu 

Consequently  It  appears  that  by  using  the  conventional  apprach  outlined 
above  one  neglects  a node  of  motion  which  may  In  fact  not  be  negligible 


Section  3. 


EXPRESSIONS  FOR  a“,  b“,  and  c“  IN  TERMS  OF  SPHERICAL 
BASISVECTORS  BUT  WITHOUT  DERIVATIVES 


We  recall  the  defintions: 


A?(k,  H)  ■ '“I  ( (U>  i 

1 U 


V' 


in  | \ ■ vj ^ f <\ 

«J<M>  - i ; f t,  **  v . * > 

\ ^ ” c)Qy  ^ ~ ; 

v A-a  v~  O l|>  ^ 


; 4 '»  -r 


C?Cy, 


V " " - ] / 't-V  - € . i £>  v1'  - 


e> 


} 


For  computation il  purposes  it  is  desirable  express  B®  and  c”  without 
derivatives  and  without  trigonometric  expressions.  We  thus  need  to 
express 


_L  oVv 

SlVv,  (?  '£)  ^ 


and 


in  such  a manner 


that  they  contain  no  derivatives.  From  the  quantum  mechanical  theory 
of  angular  momentum  (1)  one  has  the  two  operators 

V X <fc?  ^ ^ -V  v £ 0*1  & cl 


l . J - 


H & 


Vv 


l,  ^ fr~  \w  © 


These  operators  when  they  operate  on  spherical  harmonics  give 

^ \ ^ 1-  • \VA  \{  V- 4 4 \ )XJ <.  ' 


L. 


^ J V . <S  \(  U - S * ' V'1  'c  \ ' ’ 

L 


now 

and  thus 

tv" 


i e V.  - 


u 


>-Be  1 ‘v.NlK- 
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M^III  II  ijikiPWP^^jP  JIIIPBPWWWI 


or  ' rs'1' 

■ — 1 ’ 1 :.  - «.  ', 


lL-4  I ?.  |-r. 

^uhh,vumx)[-c-'v 


Cwfc 

•'D'--  ^^Xv-h^v^ve:4-\'"'-' 

so  finally  L +"  l 

ra^  - ' 


ft-  W -t  I 

\A  I 


JUlVo)  V“<_  . 


1 (lswi\) 


» v'  » 


■ 4 v-t 


LtU 


1 


- --,  i-v^vvrvi 


For  a purely  zonal  motion  one  then  has  terms  of  the  form 


C 


£ 


: LV,  t-  ;H'  . 


^ u 


V ,*  1 -s  ' U ( 

V v e J 


one 


can  note^that  this  formulation  has  no  problem  with  the  direction  at  the 

pole  since  Y1  and  Yx  are  zero  there.  The  b component  vanishes  since  it 
is  proportionate  to  t V ^ f 

c 'A 


For  purely  meridional  motion  one  has  terms  of  the  form 


U "L  b ’ L 


V ~ >Ai  * ' 'V  t / 

— s r > 


L 


1 
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Section  4.  DIRECTION  OF  THE  VECTOR  SPHERICAL  HARMONICS  AT  THE  POLES 
The  three  sets  of  vector  spherical  harmonics  can  be  written: 


^ J 6,  <f\  - e v V'u  Co.cf  ) 

(1) 

\ — ___  u 

(2) 

V _ N<v" 

C ( <>  , 4 , = - y \ 

(3) 

\)U  *-■'») 

(4) 

One  can  immediately  see  that  the  A°  have  a well  defined  direction  every- 
where on  the  sphere  i.e.,  in  the  radial  direction. 


Since  the  s and  the  C^'s  are  related  to  each  other  by 


c. 


(5) 


it  is  clear  that  M and  L being  the  same  for  each  set,  one  set  is  perpen- 
dicular to  the  other.  We  will  thus  investigate  the  behavior  of  c“  at 
the  poles.  The  behavior  of  b“'s  then  follow  form  (5). 

C™  can  be  written 


(o C\.- vu i Vi ^ 


L ' 


I' 


_ \ i l*t  h\(  L-  H * \ V 
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The  only  spherical  harmonics  which  are  non-zero  at  the  pole  are  Y°.  Con- 
sequently the  only  C®  that  can  be  non-zero  at  the  poles  are  those 


which  contain  Y^,  i.e. 


C?  * T Vf*-  'v.  -a 


_ y u up  ; ~u  o * 
- Lu|^ii|  \-e- 


V Z_  ' ' 

) V 


,2  _ „1 


Since  * 0 at  the  pole  we  only  consider  the  non-vanishing  terms, 


which  contain  Y, 


The  velocity  field  can  then  be  written  in  the  form  (at  the  pole) 

V^|Pc\o,o)Pc'  c\o,») 

• ' * 

Now,  \l  is  real 


\ » 


\f  : CL  C ^ <t  C C’ 

U L \ 


* J '+ 


4 V- 


We  have  the  following  relationship 


u- 


or 


and 


\ V 


L 


i »■ 


i l 


'Nii-L.""!*  1 _ V1 

T t - ^ 


t c c;  - 

L-  L 


or 


c c - 1 c.  . c ■ 2 c" , C . 


u v-  L.  ^ L 


\ * 

I U 


•1  ii.  - ' 

C C 

V ^ 


Thus 

• > 

<•*  \ ^ 

L C <L  * 

2 c*  c' 

Write 

-i  _ 

> U L 

V . -V  v (i,  , 

L u L 

Thus 

V 

L j U. 

— ■*  I 


*+  * ^Ll0V'|Oli 


j~  0.  L*  V " v 

-i  L~  ff  >J 


This  is  the  value  of  the  velocity  at  the  north  pole.  The  significance  is 
that  it  is  not  necessary  to  take  any  particular  computational  or  trans- 
formational steps  to  avoid  difficulties  at  the  poles  when  dealing  with 
the  velocity  field  expressed  in  terms  of  vector  spherical  harmonics. 


32 


- 


Section  5.  THE  MODEL 


This  section  will  describe  the  model  currently  implemented.  It  is  a 
global,  dry-air,  hydrostatic  balance  set  of  equations  in  spherical 
coordinates,  without  orography.  The  independent  variables  and  domains 
are: 


(the  radial  variable 


) 


\>,T  >,  O 


^ r ^ ©o 


0 J CO  latitude  O £ 6 ^ TT 

<p  j longitude  O £ <p  £ ZTT 

The  dependent  variables  are: 
v?  (<r,  e,  <p,t ) , the  horizontal  velocity  field 
, temperature  field 
, log  surface  pressure 
, geopotential 


Other  symbols  are: 


p<v,e,4>,t) . 

pressure 

Tr(r,e,<f,t)  . 

surface  pressure 

J 

acceleration  due  to  gravity 

(constant) 

-A. 

earth's  angular  velocity 

R 

gas  constant 

K 

ratio  of  specific  heats 

/N 

6-<r 

radial  unit  vector 

F(<r,e;<p,t)  , 

friction  field 

heat  forcing  function 

specific  heat  at  constant  pressure 

*“(«■/ 0/ ‘Pit)  , 

distance  from  earths  center 

radius  of  earth 

v 

\ *e  ^«Py 

&o(t)  » 

latitude  of  earth-sun  line 

longitude  of  earth-sun  line 

*p 

friction  constant 

Or) 


» heat  term  constants,  one  per  (Player. 


The  equations  are: 
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The  quantity  in  brakets  is  constrained  to  be  ^ 0. 
This  is  attempting  to  indirectly  enforce 
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The  required  initial  conditions  are: 


v<r(e,«.,0 

T (<r, 

y ( e,  <f  ,^o) 


The  boundary  conditions  are: 

$ ( 1 , 0,  <p,  t)  = O 

(1,  ^ / <f;,  t)  = f(o;0y  cf,  t)  ^ O 


The  terms  j are  in  parenthesis  to  highlight  the  fact  that 
these  quantities  are  computed  and  used  most  naturally  in  the  form.  See 
section  on  radial  functions  for  additional  discussion  of  (j/ff*)  . 
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Section  6.  PROGRAM  SYSTEM  DESIGN  CONSIDERATIONS 
INTRODUCTION 

A program  system  has  been  designed  and  implemented  to  support  research 
into  the  application  of  Vector  Spherical  Harmonics  (VSH)  to  global 
weather  models  expressed  as  sets  of  partial  differential  equations  in 
spherical  coordinates.  This  section  presents  considerations  affecting 
the  design  philosophy  and  the  tools  and  techniques  used  in  attempting  to 
carry  it  out. 

Because  the  system  is  intended  as  a general  framework  to  support  research, 
it  is  desirable  that  it  be  easy  to  use  and  easy  to  change.  A user  should 
be  required  to  specify  only  the  minimum  amount  of  information  necessary 
to  define  a given  run.  He  should  be  able  to  specify  both  simple  and 
fundamental  alternatives  - from  integration  time  step  through  integration 
algorithm  to  the  set  of  equations  to  be  solved. 

Furthermore,  using  appropriate  conventions,  many  aspects  of  code  changes 
become  mechanical.  This  is  particularly  true  of  those  changes  - storage 
structures,  loop  control  - that  distribute  over  the  system  as  the 
result  of  a change  in  a specific  part  of  the  system.  These  distributed 
changes,  if  possible,  should  be  automated  to  save  the  user  time  and 
eliminate  the  possibility  of  error.  These  considerations  imply  that 
the  system  be  highly  parameterized.  Further,  distinction  should  be 
made  between  Independent  and  dependent  parameters.  Only  those  parameters 
that  were  not  computable  from  other  parameters  would  be  independent  and 
require  user  specification. 

The  dynamic  systems  to  be  investigated  generally  require  large  amounts 
of  data  and  computation.  One  of  the  objectives  of  the  program  system 
is  to  demonstrate  VSH  advantages  in  terras  of  computational  efficiency. 

Thus  the  system,  while  being  general,  must  execute  efficiently.  Also, 
occasions  will  arise  when  efficiencies  must  be  compared  against  existing 


models  and  Implementations.  These  other  implementations  are  generally 
written  in  FORTRAN.  In  order  to  allow  the  most  direct  comparison 
of  efficiency,  and  for  other  reasons,  it  is  desired  that  the  portions 
of  the  system  applying  VSH  to  solve  equations  be  written  in  FORTRAN. 

This  gives  rise  to  a conflict.  The  high  degree  of  parameterization 
desired  for  generality,  ease  of  use  and  ease  of  change  is  not  compatible 
with  the  efficiently  executing  streamlined  code  necessary  to  best 
demonstrate  efficiency. 

The  conflict  is  resolved  using  the  PL/I  precompiler  extended  to  FORTRAN. 
Using  the  precompiler  the  system  can  be  highly  parameterized.  The 
parameters  will  be  precompiler  variables  and  the  parameterization  will  be 
resolved  at  compilation.  The  precompiler  "IF"  clause  and  "INCLUDE" 
features  allow  source  code  segments  to  be  brought  in  line  for  compilation 
by  appropriate  precompiler  variable  settings.  The  precompiler  "macro" 
feature  allows  "code  generators"  to  be  written  that  produce  FORTRAN 
code  as  a function  of  precompiler  variables.  The  precompiler  is  also 
used  to  generate  the  Job  Control  Language  procedure  to  execute  a given 
configuration  of  the  system,  insuring  consistency  between  data  definition 
statements  and  executing  code. 

VSH  Functional  Flow 


The  figure  on  the  next  page  shows  the  basic  functional  flow  of  VSH  incor- 
porating transform  methods.  Consider  the  boxes  to  be  data  resevoirs, 
or  data  sets,  and  the  numbered  lines  to  indicate  data  flow  paths  and 
a sequence  of  algorithms. 

Computations  are  done  both  in  the  physical  domain,  using  "physical" 
variables-velocity,  density,  etc-arrayed  on  layered  spherical  grids  and 
in  the  spectral  domain,  using  spectral  expansion  coefficients  of  the 
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Flow  Diagram 


physical  variables.  The  data  is  transformed  as  appropriate  to  take  advan- 
tage of  computational  and  formulational  efficiencies  in  spectral  or 
physical  space.  , 


The  flow  begins  and  ends  with  spectral  expansion  coefficients,  represented 
in  box  A. 

The  sequence  of  functional  steps  is: 

1.  Transform  coefficients  to  obtain  corresponding  values  over  the 
physical  grid  - Function  line  1 from  box  A to  box  B. 

2.  Numerically  expand  each  coefficient  over  all  layers  into  radial 
derivatives  - Function  line  2 from  box  A to  box  C. 

3.  Form  the  appropriate  spatial  derivatives  term  by  term.  - 
Function  lines  3 from  boxes  A and  C to  box  D. 

4.  Transform  to  get  values  of  spatial  derivatives  over  the  physical 
grid.  - Function  line  4 from  box  D to  box  E. 

5.  Form  the  model  equations  and  apply  boundary  conditions.  Function 
lines  5 from  boxes  b and  E to  box  F. 

6.  Transform  the  temporal  derivatives  to  coefficient  space. 

Function  line  6 from  box  F to  box  G. 

7.  Integrate  to  update  coefficients.  Function  line  7 from  box 
G to  box  A. 
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Some  of  the  areas  subject  to  change  in  the  course  of  research  include: 

1.  The  model-box  F and  function  lines  5. 

2.  Various  algorithms  - obtaining  the  initial  set  of  coefficients, 
mechanizing  the  transformations,  integration  techniques,  etc. 

3.  The  method  - the  radial  fitting  functions,  the  number  of  terms 
in  the  spectral  expansions. 

All  of  the  possible  changes  and  categories  of  change  cannot  be  anti- 
cipated apriori.  In  general,  the  primary  change  must  be  implemented 
by  recoding.  Therefore,  the  individual  algorithms  and  the  model 
equations  have  been  isolated. 

However,  a change  that  can  "analytically"  be  thought  of  as  isolated 
often  necessitates  associated  changes  distributed  throughout  the 
system. 

For  example,  adding  curl  of  curl  to  modify  the  friction  term  in  the  model 
requires  a coding  change  along  function  lines  5.  But  then  the  curl 
of  curl  must  also  be  formed  in  spectral  space  (function  line  2 and 
box  C) , and  curl  of  curl  must  be  transformed  to  physical  space  (function 
line  4 and  box  E) . The  addition  of  curl  of  curl  in  the  model  is  the 
primary,  or  independent  change  and  can  be  made  in  isolation.  The  other 
changes  follow  from,  and  are  dependent  on,  the  primary  change. 

As  another  example,  consider  changing  the  number  of  terms  in  the  spectral 
expansion  for  given  variables.  This  can  be  effected  by  setting  limits 
of  computational  loops  in  the  various  transforms  (lines  1,  4,  and  6). 

But  storage  areas  - both  core  and  direct  access  - are  also  affected. 

Note  also  that  the  first  example  implies  changes  to  various  storage  area. 
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In  both  of  these  examples  the  primary  change  requires  action  by  the  person 
making  the  change  and  is  visible  and  manageable.  The  secondary,  dependent 
changes  are  often  not  so  visible,  and  consistency  between  the  primary 
and  secondary  changes,  which  is  essentially  mechanical,  is  often  difficult 
to  implement.  An  objective  of  the  program  design  is  to  automate  where- 
ever  possible  these  secondary  changes,  relieving  the  user/changer  of  a 
major  burdern  and  assuring  consistency  throughout  the  system. 

In  short,  the  functional  flow  given  in  Figure  1 is  to  be  implemented 
to  execute  efficiently,  but  in  a manner  that  allows  for  easy  recoding 
of  primary  chanees  and  automates  the  mechanical  aspects  of  a secondary 
change. 


This  gives  rise  to  a conflict,  as  summarized  in  the  figure  below. 
Generality,  ease  of  use  and  ease  of  change  imply  a high  degree  of 
parameterization  in  the  system.  Execution  time  efficiency  implies 
minimizing  execution  time  decision  and  writing  streamlined,  highly  linear 
code.  The  implications  of  the  objectives  are  not,  on  the  face  of  it, 
compatible. 

Conflict 


Generality 

Ease  Of  Change 

Ease  Of  Use 
v i 

High  Degree  Of 
Parameterization 


Execution  Time 

Efficiency 

v / 

No  Execution  Time 
Decisions 
Streamlined  Code 


Use  of  a Precompiler 

The  solution  chosen  is  to  pay  the  price  of  flexibility  at  compile  time. 
Each  change  will  involve  a recompilation;  not  only  of  the  routine 
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containing  the  change,  but  of  all  impacted  routines.  All  control  decisions 
will  be  resolved  and  all  storage  allocations  will  be  established  at  compile 
time . 

The  tools  to  implement  this  approach  are  provided  by  the  PL/ I precompiler 
(or  preprocessor)  and  by  a prototype  FORTRAN  preprocessor  developed 
as  an  extension  of  the  PL/I  preprocessor.1  The  capabilities  of  the 
PL/l  FORTRAN  preprocessor  are: 

o Program  Libraries 

o 'Includes' 

o Macros 

o Preprocessor  Variables 

To  maximize  the  utility  of  these  tools,  all  programs  are  organized  into 
libraries  on  direct  access  storage  devices.  Libraries  are  maintained 
for  the  following  program  elements: 

o Source  Code 

o Object  Modules 

o Load  Modules 

o Job  Control  Language  Procedures 

o Linkage  Editor  Code 


1 


Developed  by  Juan  Rivero,  Wesley  Gray,  Barry  Beals  of  IBM  Federal 
Systems  Division. 


Each  library  is  organized  as  a partitioned  data  set  containing  individually 
accessible  elements,  called  members.  For  example,  an  integration  algorithm 
would  be  a member  of  the  source  library,  while  all  of  the  source  members 
necessary  to  implement  the  VSH  functional  flow  are  compiled  and  linked  to- 
gether to  produce  one  member  of  the  load  library. 

Utility  programs  have  been  developed  for  performing  housekeeping  in  the 
libraries-editing,  deleting,  adding  selecting  members. 

The  INCLUDE  capability  is  a preprocessor  feature  that  allows  source 
library  members  to  be  brought  in  line  at  compile  time.  For  example, 
the  statement  "INCLUDE  INTEG"  would  be  replaced  by  the  source  code 
contained  in  member  INTEG.  This  code  may  contain  the  statement  "INCLUDE 
DERIV"  and  similar  nesting  can  occur  indefinitely. 

The  Primary  advantage  of  the  INCLUDE  feature  is  that  the  source  code 
can  be  broken  into  logical  or  functional  units  without  unnecessary 
subroutine  linkages. 

Macros  are  PL/ I procedures  that  execute  at  compile  time.  Their  chief 
attraction  is  as  code  generators.  A macro  can  return  FORTRAN  code 
that  is  functionally  dependent  on  the  input  argument  values.  Macros 
combined  with  preprocessor  variables  are  the  means  by  which  distributed 
changes  are  automated. 

Preprocessor  variables  are  PL/I  variable  that  are  active  at  compile 
time.  They  are  of  two  types,  character  strings  and  integers.  They  are 
used  in  preporcessor  and  non-preporcessor  statements.  Preprocessor 
statements  (assignment,  if  then  else,  do  loops)  are  executed  at  compile 
time.  Preprocessor  variables  are  used  in  these  statements  in  the  same 
way  that  PL/I  variables  are  used  in  PL/I  statements.  The  Preprocessor 
IF  statement  allows  selective  execution  of  INCLUDE  statements. 


Where  preprocessor  variables  appear  in  non-preprocessor  statements 
(FORTRAN  code)  they  are  replaced  by  their  current  value. 

For  example,  the  dimension  of  a FORTRAN  array  may  be  coded  with  a 
preprocessor  variable  in  place  of  a literal  size  specification.  When 
the  preprocessor  encounters  this  statement  the  varible  will  be  replaced 
with  its  current  value.  The  following  figure  gives  an  example  of  a pre- 
processor macro  that  produces  FORTRAN  code. 

MACRO  CALL 

^EQUIVALENCE  (V//B//  ,2,3,  $DIMENS IONS  s WORK , 201) 

GENERATED  CODE 

EQUIVALENCE (V1B1(1) ,W0RK(201)) , (V1B2(1) ,W0RK(301)) , 

X (V1B3(1) , W0RK(401) ) 

EQUIVALENCE  (V2B1(1) ,WORK(201)) , (V2B2(1) ,WORK(3201) ) , 

X (V2B3(1),W0RK(6201)) 

@Equivalence  is  the  macro  name.  The  arguments  V//B//  and  WORK  are  literals. 
The  remaining  four  arguments  can  be  literals  or  preprocessor  variables. 

In  actual  use  they  are  usually  preprocessor  variables  whose  values 
have  been  set  from  some  compile  time  control  parameters. 

Storage  is  arranged  in  blocks  according  to  the  function  being  performed. 
Naming  conventions  have  been  established  for  the  blocks.  In  this  example, 
V//B//  is  the  "name  base"  for  allocating  blocks  of  storage  for  vector 
variables.  The  second  and  third  arguments  in  the  macro  call  give  the 
number  of  vector  variables  and  the  number  of  storage  blocks  to  be  allocated 
to  each  variable.  There  are  thus  two  vector  variables  and  each  has  three 
blocks.  The  $DIMENSI0NS  string  has  two  entries  - 100  and  3000.  The  size 
of  all  blocks  for  vector  variable  1 will  be  100  words  and  the  size  of 
all  blocks  for  vector  variable  2 will  be  3000  words. 
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Allocation  and  alignment  of  storage  - the  number,  size,  and  relative 
locations  of  various  blocks  - can  be  automated  using  similar  naming 
conventions  for  other  types  of  storage  blocks.  Also,  the  preprocessor 
variables  that  act  as  parameters  to  the  macro  calls  can  themselves 
be  computed  at  compile  time  from  some  minimal  set  of  independent  pre- 
processor variables.  This  double  level  of  automation  reduces  the  load 
on  the  user /changer  and  assures  consistency. 


Section  7.  PROGRAM  SYSTEM  FUNCTIONAL  DESCRIPTION 


The  main  function  of  the  program  system  is  to  apply  Vector  Spherical 
Harmonics  (VSH)  to  global  weather  modeling.  This  is  done  in  the  PREDICT 
module.  This  implies  the  need  for  several  support  functions  such  as 
initialization,  (the  Input  module)  report  generation  (the  Output  and 
Polar  modules)  and  tabular  data  construction  (the  TABLE  module). 

The  weather  model  currently  implemented  is  a Mintz-Arakawa  dry  air 
formulation  in  sigma  coordinates  on  eight  layers.  The  Richardson  two 
step  scheme  is  used  for  time  integration. 

Since  the  major  portion  of  the  computations  are  done  with  spectral 
expansion  coefficients,  an  interface  - the  Standard  Format  - has  been 
developed.  Initialization  data  derived  from  any  source  must  be  put  into 
the  Standard  Format  before  being  passed  to  the  PREDICT  module.  Similarly, 
the  time  history  generated  by  PREDICT  is  in  the  Standard  Format,  which 
must  be  converted  into  desired  report  format. 

The  current  source  of  initialization  data  is  the  National  Meteorotogical 
Center  (NMC)  GLOPEP  data  set.  This  is  "live"  global  meteorological  data 
updated  on  a daily  basis,  and  thus  it  also  provides  a means  of  evaluating 
forecasts  produced  by  PREDICT. 

System  Components 

The  program  system  consists  of  the  following  major  components: 


o Input 

Converts  data  from  GLOPEP  format  to  an  initial  Mintz-Arakawa 
model  state.  Also,  on  option,  selects  an  initializing  state, 
or  set  of  states,  from  a Standard  Format  PREDICT  time  history. 
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o 


Predict 


Applies  VSH  techniques  to  a Mintz-Arakawa  cast  of  a global 
weather  model  to  produce  predictive  estimates  of  weather  be- 
havior . 


Converts  selectM  samples  from  a PREDICT  time  history  in 
Standard  Format  into  GLOPEP  format  for  subsequent  mapping  and 
analysis.  Prints  user  selected  data  in  various  spectral 
and  physical  representations  from  the  PREDICT  Standard 
Format  time  history. 


o Polar 


Produces  user  selected  polar  stereographic  contour  plots  of 
GLOPEP  data. 

o Table 

Prepares  various  time  invariant  data  tables  used  in  the 
VSH  protions  of  PREDICT. 

Input 

An  initial  state,  or  states,  must  be  prepared  for  PREDICT.  There  are 
two  possible  sources  - GLOPEP  or  a PREDICT  time  history.  The  state 
must  contain  the  following  data  in  spectral  expansion  coefficient  form. 

o wind  velocity  ( 0" ) , on  all  layers 

o temperature  over  sigma  (//(T)  on  all  layers 

o natural  logarithm  of  sea  surface  pressure  (^) 

48 


If  the  source  of  the  initializing  data  is  a PREDICT  time  history,  the 
necessary  variables  are  simply  selected  at,  or  within  a specified  tolerance 
of,  the  desired  time.  If  the  source  of  the  initializing  data  is  a GLOPEP 
data  set,  conversions  must  be  made  to  obtain  the  variables  Tir  and 


This  is  a three  step  process: 

o Convert  GLOPEP  variables  to  PREDICT  variables 

o Convert  from  GLOPEP  radial  layers  to  PREDICT  radial  layers 

o Convert  from  physical  variables  on  latitude  longitude  grids 
to  spectral  expansion  coefficients. 

Applicable  GLOPEP  data  consists  of: 

U = horizontal  component  of  wind  from  West  to  East 

V = horizontal  component  of  wind  from  South  to  North 

T = temperature 

Z = pressure  height 

This  data  is  given  in  pressure  coordinates  on  concentric  spherical  grids, 
equally  spaced  in  latitude  (©)  and  longitude  (<$>)  at  2.5°.  The  grids  are 
set  at  pressure  heights  of  1000,  850,  700,  500,  400,  300,  250,  200,  150, 
100,  70,  50  millibars. 
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o 


Conversion  to  PREDICT  variables 


GLOPEP  velocity  is  transformed  to  rectangular  coordinates  at  each 
GLOPEP  (&$)  point  by 
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GLOPEP  temperature  is  simply  copied  at  this  step.  It  will  be  divided 
by  sigma  after  conversion  to  spectral  expansion  coefficients. 

Sea  surface  temperature  is  currently  taken  to  be  1000  mb  temperature. 
It  may  optionally  be  taken  as  GLOPEP  surface  temperature,  but  this 
includes  orography,  and  the  current  Mintz-Arakawa  model  is  a spherical 
model. 


Sea  surface  pressure,  /C  , is  obtained 


condition 


using  the  hydrostatic  balance 


where  p is  pressure,  r is  radius,  g is  gravity,  T is  temperature  and 
R is  the  universal  gas  constant. 


Since  sea  surface  is  very  near  1000  mb  pressure,  one  takes  the  slope 
of  pressure  versus  r at  1000  mb  from  the  hydrostatic  balance  condition 
to  use  in  a point  slope  form  for  linear  extrapolation  from  1000  mb 
to  sea  surface  pressure 

7C  - Ictoo  ^ ~ 

where  re-rlOOOmb  = Z lOOOmb,  the  first  layer  GLOPEP  pressure  height. 
Adjusting  for  conversion  from  GLOPEP  to  PREDICT  units, 

7c  4 lo1 

and 

P - 
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Conversion  to  PREDICT  radial  layers. 


The  GLOPEP  data  are  given  at  the  12  pressure  heights  listed 
previously.  PREDICT  requires  data  on  a number  (currently 
8)  of  Gaussian  spaced  layers  corresponding  to  the  number  of 
orthogonal  polynomials  used  in  radial  expansions.  Data  for 
the  PREDICT  layers  are  derived  from  the  given  data  nn  the 
GLOPEP  layers  by: 


Define  a function  on  the  GLOPEP  data  by  fitting  straight 
lines  between  the  GLOPEP  data  points  . 


o Expand  this  function  in  orthogonal  polynomials 

o Evaluate  this  expansion  on  the  PREDICT  layers, 

o Conversion  from  physical  to  spectral  representation. 


As  detailed  elsewhere  in  this  report,  spectral  expansion 
coefficients  xml  for  given  m,  1 truncation  limits  are  obtained 
from  physical  data  x(©;<$>)  by 
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Integration  over  <P  is  by  the  trapezoidal  rule,  while  integration  over 
O'  is  by  Simpson's  rule  with  end  point  correction. 


PREDICT 


The  PREDICT  module  generates  numerical  solutions  of  a global  weather  model. 
The  current  weather  model  is  a Mintz-Arakawa  hydrostatic  balance,  dry  air, 
smooth  earth  model  expressed  in  sigma  coordinates.  It  is  described  in 
greater  detail  elsewhere  in  this  report. 


PREDICT  has  been  built  to  apply  VSH  techniques  to  the  solution  of  weather 
models.  These  techniques  have  been  described  in  general  in  "Weather  Systems 
Department  IRAD  Report  1972"  and  in  application  to  the  Mintz-Arakawa  model 
in  "Global  Weather  Modeling  with  Vector  Spherical  Harmonics  Report  No.  2" 
ARPA  order  number  2089,  program  code  number  G270GE. 

Brief  descriptions  of  the  PREDICT  functional  components  will  be  given  here, 
along  with  a somewhat  more  detailed  exposition  of  the  steps  taken  to 
form  the  Mintz-Arakawa  model. 


The  organization  of  PREDICT' s functional  components  is  shown  in  the 
diagram  below. 
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Initialization 


The  initialization  section  reads  in  the  initial  state,  or  states,  for  the 
model,  and  various  constants  and  tables  and  run  control  data,  such  as 
integration  time  step  and  time  span.  Requisite  data  not  supplied  as 
input  are,  where  possible,  given  default  values.  After  various  validity 
and  consistency  checks,  control  is  passed  to  the  time  integration  module. 

Time  Integration 

PREDICT  currently  uses  Richardson's  method 

for  time  integration.  This  is  a two  step  method.  When  only  one  initial 
state  has  been  input,  a standard  fourth  order  Runge-Kutta  integration 
is  used  to  provide  the  other  starting  value.  The  program  is  written 
so  that  Runge-Kutta  can  optionally  be  used  as  the  only  integration  scheme, 
but  this  option  has  not  to  date  been  exercised  due  to  the  computational 
cost  of  Runge-Kutta. 

All  time  integrations  are  performed  in  spectral  space.  For  each  time  up- 
date, the  integration  section  calls  the  model  computation  routine  to  form 
the  temporal  differential  equations. 

Model  Computations 

The  sequence  of  steps  to  form  the  Mintz-Arakawa  model  is  outlined 
below.  This  sequence  is  performed  four  times  for  a Runge-Kutta  inte- 
gration step  and  once  for  a Richardson  integration  step. 

In  spectral  space: 

o Do  radial  truncation  on  U~,TjS7  This  involves  fitting  data 
computed  by  the  previous  integration  step  with  the  number  of 
polynomials  used  to  represent  radial  behavior.  It  is  done 
to  avoid  aliasing  in  subsequent  computations. 
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o Set  radial  component  of  velocity  to  zero  on  all  layers,  assuming 
only  horizontal  velocities. 

o FormV^on  layer  1 and  transform  to  physical  space. 

On  all  layers: 


o Compute  the  components  of  kinetic  energy  for  this  layer  and 
sum  into  total  kinetic  energy  for  this  time. 

o Transform  Lr  to  physical  space. 


o Compute  ^7*  Cr-  and  (tr.v)(Tand  transform  to  physical  space. 

o If  not  layer  1,  where  (l/(T)  is  held  constant,  form  V (~^/^) 

and  transform  it  and  (T/Vr)  to  physical  space. 


Form  ^*/^Oin  spectral  space. 


o Form  (y~-V^  and  the  integral  equation  for 


in  physical  space. 


Transform  U~.  V^from  physical  space  to  spectral  space  and  back 
to  physical  space  for  proper  truncation  to  avoid  aliasing 
when  forming  the  triple  product V^in  the  JCT/rfrt  equation. 
This  is  done  since  the  physical  grid  is  sized  for  simple  products. 


o Transform  the  integral  for 


to  spectral  space. 


In  spectral  space: 


o Compute  0 on  all  layers  by  radial  quadrature. 

o Compute  on  layer  1 by  radial  quadrature  and  transform 

to  physical  space. 


o 


Form  on  all  layers,  solve  for  by  radial  quadrature 

and  transform  G~/(p  to  physical  space. 


Compute  ^ /u6~' and 
physical  space. 


on  all  layers  and  transform  to 


o Compute  S7<?  on  all  layers  and  transform  to  physical  space. 


In  physical  space: 


Form 


on  all  layers  and  transform  to  spectral  space. 


o Form  on  all  layers  but  the  first  and  transform  to 

spectral  space. 


In  order  to  perform  these  computations,  the  model  calls  on  the 
Differential  Operators,  Radial  Quadrature  and  Transforms  sections. 


Differential  Operators 


This  program  section  performs  vector  and  scalar  differential  operations 
on  given  variables  according  to  formulae  given  by  VSH.  For  the  Mintz- 
Arakawa  model,  only  horizontal  differential  operators  are  used,  so  no 
radial  derivatives  appear.  The  various  operators-curl , divergence, 
gradient,  advection  - are  formed  as  combinations  of  the  numerical  con- 
stants e^,^,  l/r  and  the  spectral  expansion  coefficients.  The  exact 


VSH  formulation  for  those  operators  is  given  in  "Global  Weather  Modeling 
with  Vector  Spherical  Harmonics,  Report  No.  2"  referenced  earlier  at  the 
beginning  of  the  discussion  of  the  PREDICT  module. 


In  the  event  of  a three  dimensional  model,  radial  derivatives  of  the 
expansion  coefficients  w">uld  be  used  in  forming  the  differential  operators. 
These  would  be  supplied  by  the  Radial  Expansion  section. 
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Transforms 


This  section  performs  the  transformations  of  quantities  from  spectral 
expansion  coefficients  to  physical  grid  representations  and  vice-versa. 
For  computational  efficiency,  the  model  forms  equations  involving 
products  in  physical  space,  first  computing  the  individual  terms  in 
spectral  space  and  then  transforming  them.  The  resulting  equation 
is  then  transformed  back  to  a spectral  representation  for  a subsequent 
integration  or  quadrature. 


The  form  for  transforming  from  physical  to  spectral  space  is 

r*r** 

--  J I x<>,  <p)  yjr-  %,  •, ; a s- 


O o 

The  form  for  transforming  from  spectral  to  physical  space  is  given 
by  the  truncated  harmonic  expansion 


M -I 


K( 


L-l 

y i ^ 


The  transforms  and  the  exact  techniques  used  are  discussed  in  greater 
detail  elsewhere  in  this  report. 

Radial  Quadrature 

The  solution  of  the  Mintz-Arakawa  model  involves  integrations  over 

the  radial  variable  sigma.  This  is  done  by  Gaussian  quadrature  on  spectral 

expansion  coefficients.  The  PREDICT  radial  layers  are  Gaussian  spaced. 

The  initial  data  prepared  in  Input  is  derived  radially  using  the  orthogonal 
polynomial  set  to  be  used  in  PREDICT.  This  assures  that  the  PREDICT 
data  is  exactly  representable  radially  by  a polynomial  of  highest  degree 
consistent  with  the  number  of  layers  in  PREDICT. 
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Radial  Expansion 


J(ZS  c'^/u  /A, 

In  the  Mintz-Arakawa  model,  the  equations  for  /£xand  /j'L.  require 
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the  terms  j^-and  The  radial  derivatives  of  (r  and  ~f~/(T  are 

computed  in  spectral  space  by  analytically  differentiating  the  radial 
function  expansions  of  the  variables.  This  is  possible  since,  as  mentioned 
above,  the  data  are  derived  to  be  exactly  representable  by  the  orthogonal 
polynomial  set,  currently  Legendre  Polynomials  of  the  first  kind,  which  are 
differentiable. 


Predict  Output 

At  the  conclusion  of  the  time  step,  the  Integration  section  calls  the 
PREDICT  OUTPUT  section  (not  to  be  confused  with  the  OUTPUT  program, 
which  is  a separate  entity).  Here  various  user  selected  data  are 
printed,  and  the  Standard  Format  time  history  data  set  containing  user 
selected  data  at  user  selected  intervals  is  prepared. 

Output 

The  Output  module  has  two  main  functions: 

o Print  selected  data  in  selected  formats  from  the  Standard 
Format  time  history. 

o Prepare  a GLOPEP  data  set  at  a selected  time  instant  from  data 
on  the  Standard  Format  time  history. 

Print 


Any  data  on  the  time  history  data  set  may  be  printed.  Fields  may  be 
printed  in  the  following  forms  an  selected  layers 
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Scalar  Spectral: 


Vector  Spectral 


Vector  Physical: 


Scalar  Physical 


Coefficients  of  Y® 

Coefficients  of  A®,  B®,  C® 
Tm 

jl 

at  grid  points  in 

components 
components 
X,Y,Z  components 

at  grid  points 


No  capability  currently  exists  in  OUTPUT  to  transform  fields  from  physical 
to  spectral  space  or  vice  versa  for  printing.  The  only  transforms  done 
are  from  one  representation  to  another  in  the  given  space.  The  Standard 
Format  time  history  carries  fields  in  the  following  forms: 


Scalar  Spectral: 
Vector  Spectra]: 
Vector  physi;al: 

Scalar  physical: 

Glopep 

Output  can  prepare  a GLOPEP  data 
in  three  steps: 


Coefficients  of  Y® 

Coefficients  of  T®^ 

X,Y,Z  coordinates  atf<9;(p) 
grid  points 

at  (&,$)  grid  points 

set  at  a selected  time.  This  is  done 
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o Transform  PREDICT  spectral  fields  to  physical  data  on  GLOPEP 
cf  grid  at  PREDICT  layers  . 

o Transform  from  PREDICT  layers  to  GLOPEP  layers, 
o Transform  from  PREDICT  variable  to  GLOPEP  variables. 


Spectral  to  Physical  Transformation 


Physical  data  is  computed  from  spectral  expansion  coefficients  X™ 


by 


,1 


Radial  Transformation 

Obtaining  data  on  GLOPEP  layers  from  data  on  PREDICT  layers  involves  two 
steps: 

o Determine  the  coefficeints  of  a Legendre  Polynomial  expansion 
of  the  data  on  the  PREDICT  layers  by  Gaussian  Quadrature. 

o Evaluate  the  Legendre  Polynomial  expansion  on  the  Glopep 
layers . 

Transform  to  GLOPEP  Variables 

The  GLOPEP  data  set  requires  the  variables: 

U horizontal  component  of  velocity  from  West  to  East 


V 


horizontal  component  of  velocity  from  North  to  South 


Z Pressure  height 

T Temperature 

TSURF  Surface  Temperature 

PTROP  Tropopause  Pressure 

RH  Relative  Humidity 

These  data  are  obtained  from  PREDICT  data  as  follows: 
o Velocity 


u 

as 

- sin  $ cos  <P 

0 

V 

- cosmos®"  -sin<pcos 

sin 

o Pressure  height 

a - 0/g 

where  0 is  PREDICT  geopotential  and  g is  gravity 

o Temperature  is  simply  set  equal  to  PREDICT  temperature  that 

has  been  carried  thru  the  previous  two  transformations 

o Surface  temperature  is  set  equal  to  1000  millibar  temperature 

o Tropopause  pressure  is  made  to  carry  PREDICT  surface  pressure 

for  subsequent  mapping  via  POLAR.  Relative  humidity  is  set  to 
zero. 

Polar 


Polar  generates  contour  maps  of  selected  data  from  a GLOPEP  data  set. 
User  controlled  parameters  allow  control  of  such  things  as  map  size, 
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data  scaling,  contour  interval,  hemisphere,  layer  and  variable  selection. 

The  program  defaults  to  produce  hemispheric  polar  stereographic  projections, 
but  any  desired  projection  may  be  generated  by  suitable  choice  of  control 
parameters. 


Table 


TABLE  generates  various  time  invariant  data  required  by  PREDICT.  Among 
these  data  are: 


o 


o 


and  F® 


r sets  for  use  in  differntial  operators 
sets  for  use  in  transforms 


Matrices  of  Legendre  Polynomials  with  appropriate  factors  for 
use  in  radial  expansion  and  radial  guadrature 
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Section  8.  RADIAL  FUNCTION  CONSIDERATIONS 


The  first  program  runs  exhibited  disappointing  reconstructions  of  initial 
condition  data  fields.  This  was  primarily  due  to  poor  radial  expansions, 
described  here.  See  also  the  section  on  vector/scalar  spherical 
harmonic  expansions. 

Briefly  the  situation  was: 


1. 


2. 


Given  data  fields,  on  a 2.5°  latitude,  longitude  grid  on 
twelve  pressure  surfaces  (1000,  850,  700,  500,  400,  300, 
250,  200,  150,  100,  70,  50  mb.).  The  NMC  at  Suitland,  Md. 
kindly  supplied  these  data  fields. 


Expand  the  fields  in  truncated  expansions  of  Legendre  poly- 
nomials time  vector /scalar  spherical  harmonics 


P^(S)  Y*(0,  <p) 


°r  p*(s)Tire'<f) 


Note  that  s is  a linear  mapping  of  <r  (our  true  radial  variable) 
such  that  l + is  equivalent  tc  (the  interval  of 

orthogonality  of  Legendre  polynomials).  The  € is  because  the 

1 ,west  lflyer  for  Gausian  quadrature  over  s is  made  to  correspond 
to  (Ta  | . See  Figure  8-1. 


62 


'pww'r^'irr^^j.yp^f  1 '■  h rvvjs.  W!*5 w wrr»“ 


NMC  Initial 
Data  Layers 


Gaussian  Computation 
Layers  (5  Layer  Case) 


Gaussian  Computation 
Layers  (8  Layer  Case) 
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Figure  8-1.  Radial  layers. 
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3.  Evaluate  the  spectral  expansions 

( 51  /C-CA^-fc)  R,  (s)Vj * (e,  <q) 

field  = < 


f^(s)  T a ($,  «) 

at  the  grid  points  of  the  given  data  fields. 

4.  Compare  given  field  vs.  reconstructed  fields. 

These  first  program  runs  performed  radial  expansions  truncated  to  four 
Legendre  polynomials  (Pq  thru  P^.  This  requries  five  Gaussian  spaced 
layers  over  “1$  I to  handle  product  non-linearities  without  aliasing 
In  general,  expansions  in  Pk>  for  k - 0,1,  . . .K,  the  minimum  number  of 
layers  (to  avoid  aliasing  with  a product  nonlinearity)  is 


K 

Number  of  Lavers 

3 

r 

5 

4 

7 

5 

8 

6 

10 

7 

11 

8 

13 

9 

14 

10 

16 

11 

17 

The  most  isolated  major  discrepancy  in  the  reconstructed  data  fields 
involved  temperature.  In  this  case,  step  2 above  included  dividing  the 
given  temperature  field  (T)  by  <T,  then  expanding  (T/r) . Step  3 included 
a multiplication  by  <T  to  recover  the  desired  field.  The  problem  was 


64 


simply  that  (T/fl-)  could  not  be  fit  well  enough  (radially)  with  a reasonable 
number  of  polynomials,  even  though  T could  be.  The  sketch  in  Figure  8-2 
illustrates  this  problem.  The  solution  adopted  was  to  expand  the  field 
T instead  of  (T/<r) . In  implementing  the  change  it  was  easiest  to  form 
T/<r  and  use  this  quantity  in  the  computations  even  though  the  ultimate 
spectral  expansion  is  now  of  the  field  T. 

The  quantity  (T/<r)  was  originally  chosen  for  two  reasons. 


a. 


Implicit  upper  boundary  condition  on  T.  - That  is  for  any  m,l 
tne  coefficient  in  the  expansion  of  (T/<r)  is  of  the  form 


K 


This  polynomial,  when  evaluated  at  s = 1 


(corresponding  to  (T®  0,  t~  =oO)  will  have  a possibly  large  but 
finite  value.  Hence  T at  f=o  must  be  equal  to  zero.  The 
current  expansion  of  T removes  this  boundary  condition. 


b.  Maximum  aliasing  control.  - With  expansion  of  (T /<r)  the  sources 
of  aliasing  error  are  limited  to  the  forcing  functions  (friction 
and  heating  terms)  and  the  application  of  the  lower  boundary 
condition  on  temperature.  Now  the  equation  of  motion  is 
subject  to  radial  aliasing  through  the  term  V$  . 


We  then  studied  radial  variations  for  several  time  samples  of  NMC  initial 
condition  data.  This  was  done  by  expanding  the  data  (temperature, 
pressure  height,  velocity)  in  spherical  harmonics  on  the  NMC  pressure 
surfaces.  Examinations  of  the  radial  variations  of  particular  harmonic 
coefficients,  with  significant  magnitudes,  revealed  three  characteristic 
shapes.  These  are: 
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given  temperature 


T gemperature  recovered  from  the  fit  of  — 

' a 


Figure  8-2.  Temperature  fitting  problem. 


J 
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a. 


The  Y°  coefficient  of  pressure  height.  Figure  8-3. 
o 

b.  The  Y°  coefficient  of  temperature.  Figure  8-4. 

o 

c.  Most  other  coefficients,  real  & imaginary  parts,  (e.g.,  the 
Y°  coeficient  of  pressure  height.  Figure  8-5.) 

The  remaining  coefficients  tended  to  be  either  of  insignificant  amplitude 
(in  which  case  measurement  or  computational  noise  could  be  the  cause 
of  the  radial  structure  variation)  or  more  smooth  than  (c)  above.  Be- 
cause of  their  dominance,  cases  (a)  & (b)  are  singled  out  even  though 
they  have  less  extreme  variation  over  (T  then  (c) . Figures  8-3,  8-4 
and  8-5  illustrates  typical  values  of  these  coefficients  and  their 
reconstruction  after  expanding  in  P^(S)»  k = 0,1,... K for  selected 
values  of  K. 


Observe,  in  Figure  8-3,  the  relatively  good  expansion  with  K - 3, 

of  the  Y°  coefficient  of  pressure  height.  But  even  for  moaerateiy 
o 

large  K,  the  errors  are  very  significant  because  of  the  extreme  magnitude 
of  this  coefficient  (compare  with  the  next  largest  coefficient  for  this 
field,  Y°,  in  Figure  8-5).  These  errors  are: 


K 

Maximum  error 

in  Y° 

0 

<T  at  which  max 
error  occurs 

average  error 

over  12  sample  points 

3 

1.02  km 

.10 

.47  km 

5 

.37 

.10 

.14 

7 

.25 

.05 

.07 

Now  in  our  model,  pressure  height  (actually  geopotential,  $ ) reenters 

the  computations  only  via  7$  in  the  equation  of  motion.  It  so  happens 

that  Y°  has  no  horizontal  component,  hence  our  prediction  process  is 
o 

not  corrupted  by  this  sizable  error.  The  error  will  ultimately  be 
fixed  by  expanding  (for  this  coefficient  only)  the  deviation  from 
the  nominal,  computing  this  quantity  over  time,  and  then  adding  the 
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Figure  8-3.  Radial  structure  of  Yq  coefficient  of  pressure  height 
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deviations  to  the  nominal  for  output.  For  the  time  being  our  maps 
of  pressure  height  contain  this  error  which  appears  as  a bias  at 
each  output  level.  Typical  values  of  this  bias  error  (K-5)  are: 


Pressure 

Level (True  height  - computed  height) 


1000  mb 

-87 

850 

3 

700 

67 

500 

-197 

400 

-113 

300 

118 

250 

186 

200 

140 

150 

- 64 

100 

-374 

70 

-365 

50 

17 

In  Figure  8-4,  we  see  an  acceptable  fit  for  K =*  3 and  a good  fit  for 
K = 5,  with  respect  to  the  Y°  coefficient  of  temperature.  The  errors 
are: 


(T  at  which  max  average  error 

K Maximum  error error  occurs over  12  points 


3 

5.42 

°K 

.10 

2.48 

°K 

5 

2.28 

°K 

.10 

.77 

°K 

7 

1.02 

°K 

.05 

.56 

°K 

71 


ttiii 


Figure  8-5  shows  several  expansions  of  the  Y°  coefficient  of  pressure 
height.  The  independent  axis  labels  (and  errors  below)  represent  the 
situation  at  the  poles,  where  Y^  assumes  it's  greatest  value.  The  error 
summary  for  this  coefficient, 


K 

Maximum  error 

IT  at  which  max 
error  occurs 

average  error 
aver  12  points 

3 

192  m. 

.15 

83  m. 

5 

67  m. 

.07 

26  m. 

7 

51  m. 

.15 

20  m. 

It  is  interesting  to  note  that  the  (K=5)  expansion  is  superior  to  the 
(K=7)  expansion,  for  this  coefficient,  over  the  lower  layers.  The 
errors  in  meters,  by  layer  are: 

<T  K — ► 3 5 7 

i 


1.0 

32 

9 

15 

.85 

71 

9 

10 

.70 

23 

4 

16 

.50 

95 

3 

11 

.40 

91 

16 

0 

.30 

3 

21 

26 

.25 

82 

5 

12 

.20 

169 

46 

28 

.15 

192 

56 

51 

.10 

69 

18 

1 

.07 

52 

67 

42 

.05 

110 

53 

32 

Based  on  these  investigations  it  appeared  reasonable  to  start  our  prediction 
runs  using  radial  expansions  up  to  and  including  P^,  which  requires  eight 
computational  layers.  Clearly  cubic  expansions  are  not  adequate,  and  the 
improvement  of  seventh  degree  expansions  seems  small  compared  to  the  38% 
increase  of  computational  effort  over  the  fifth  degree  expansion  (i.e., 

8 layers  ->  11  layers). 

After  making  the  changes  described  here,  our  data,  reconstructed  after 
three  dimensional  expansion,  was  in  much  better  agreement  with  the  source 
fields.  Except  for  pressure  height  bias,  we  feel  that  the  bulk  of  any 
discrepancies  is  due  to  fairly  low  spherical  harmonic  truncation  limits 
we  have  used  to  date. 


Section  9.  VECTOR/ SCALAR  SPHERICAL  HARMONIC  EXPANSIONS 


In  our  spectral  method,  a scalar  field  is  approximated  by  a three- 
dimenBional  truncated  spectral  expansion, 

K N L,  ^ 

* v,  1')  = ^ ft  cs)  y,  (e,  «o  ; 


-^=0  >''  = 41-  I 

and  a vector  field  by, 
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The  radial  expansions  are  discussed  in  another  section  of  this  report. 


Here  we  will  address  only  the  horizontal  expansions 

in  l. 
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The  scalar  expansion  is  the  key  item,  since  the  vector  transformations 


are  handled  as  follows.  Represent  the  vector  field 
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Perform  a triple  scalar  transform  on  the  components 
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These  components  are  then  mapped  to  the  T^  basis  by 
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The  vector  spherical  harmonic  coefficients  are  then  obtained  using  the 
relationship 
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and  its  inverse 
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Now  let's  examine  some  of  the  details  of  a scalar  expansion.  The  expansion 
coefficients  are  found,  from  the  orthonormality  of  the  Y^'s  by 
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This  is  implemented  by  integrating 
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followed  by 
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During  the  prediction  steps,  these  quadratures  are  exact  since  the  fields 
are  truncated.  The  «p  integration  is  trapezoidal  and  the  6 integration 
is  Gaussian.  For  expanding  initial  condition  data,  supplied  on  a 
uniform  mesh  0/<p  grid,  we  first  tried  trapezoidal  quadrature  for  both 
integrations.  This  resulted  in  significant  errors  near  the  poles.  The 
errors  were  removed  by  using  Simpsons  quadrature  with  end  point  correction 


for  the  6 integration.  This  was  easy  to  do  because  of  the  sinB  factor. 
That  is,  end  point  correction  requires 
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and  no  explicit  differentiation  is  required 
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In  parallel  with  prediction  runs,  made  with  fairly  low  truncation  limits 
(M=L=9),  we  studied  the  effect  of  horizontal  truncation  limit  variation 
on  a scalar  field  (1000  mb  pressure  height). 


For  a real  field,  one  doesn't  need  to  compute  the  expansion  coefficients 
with  negative  m,  since 

c (1 ,-m)  = (-1)“  c (l,m)* 


and  for  real  vector  fields 
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Hence  we  can  represent  our  expansion  coefficients  schematically  as  the  shaded 
portion  of: 


We  have  chosen  to  implement  a truncation  policy  that  can  be  called 
trapezoidal.  That  is,  given  M,L  (L^.  M)  we  retain  the  coefficients  indicated 


\ 
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Other  spectral  modelers  seem  to  prefer  a rhombiodal  truncation  policy, 
based  on  two  parameters  M,J.  This  truncation  is  illustrated  by 


Now  let's  list  some  of  the  factors  associated  with  variations  of  the  trun- 
cation parameters  of  the  two  truncation  policies. 


a.  Number  of  coefficients 

b.  Minimum  number  of  uniformly 
spaced  points  required  to 
transform  a product  non- 
linearity without  aliasing 

c.  Minimum  number  of  Gaussian 
spaced  0 points  required  to 
transform  a product  non- 
linearity without  aliasing 


Trapezoidal 


Rhombiodal 
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d.  Approximate  operation 
count  for  the  Legendre 
step  in  a Y™  trans- 
form, either  direc-  (3L+1)  (M+l)  (2L-M+l)/4  (J+l) (M+l) (3J+2M+l)/2 

tion.  (Minimum 
physical  grid). 

The  Fourier  part  of  a Y™  transform  can  be  implemented  in  various  fashions 
with  both  truncation  policies.  Its  operation  count  could  range  from  one 
to  a small  fraction  times  (d)  above. 

Let's  now  look  at  two  ways  to  define  equivalence  between  the  truncation 
policies,  and  compare  the  'quality'  of  equivalent  truncated  expansions  of  a 
particular  field.  This  field  (1000  mb  pressure  height,  hour  0,  25  July  1974) 
was  used  only  because  it  was  the  field  with  the  most  'structure'  on  the 
first  data  tape  we  received  from  the  National  Meteorological  Center. 

Because  of  the  similarity  of  the  parameter  M,  we  shall  require  equivalent 
truncation  limits  to  have  the  same  value  for  M.  One  possible  way  to  relate 
the  parameters  L and  J would  be  to  produce  the  same  number  of  coefficients. 

That  Is  from  (a)  above, 


or 


(M+L ) ( 2L-M+2 ) / 2 
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Another  possiblity  is  to  equate  the  operation  counts  given  in  (d)  above, 
since  the  spectral  4 — ► physical  transformations  dominate  the  computation 
effort  in  spectral  methods.  In  this  case 
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For  12^.M^20  and  M<L^40 
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We  will  split  the  difference  and  say  that 
J ■ L - (M+l)/2  are  equivalent. 

One  measure  of  the  truncation  error  can  be  formed  by  the  sum  of  the 
magnitudes  squared  of  the  unretained  coefficients.  This  is  proportional 
to  the  surface  integral  of  the  error  field  squared. 

1 ,m 

unretained 

where  )^is  t*le  truncated  expansion  of  X. 

Figure  9-1  is  60  meter  contour  plot  of  the  untruncated  1000  mb  pressure 
height.  Note  that  this  plot  (and  the  others  in  this  section)  do  not  quite 
wrap-around  in  tp  . 

Figure  9-2  tabulates  the  coefficient  magnitudes  (1  and  m up  to  25)  obtained 
by  expanding  the  1000  mb  pressure  height.  These  have  been  multiplied  by 
10  and  rounded  to  integers.  Note  that  the  RSS  and  RMS  figures  are  with 
respect  to  the  rows  and  columns  of  the  triangular  array  printed  except 
that  the  M»L=0  term  has  been  excluded.  That  is,  the  factor  of  10  is  in- 
cluded and  the  implied  coefficients  (with  m negative)  have  not  been  in- 
cluded in  the  RSS,  RMS  computations. 

Figures  9-3  show  the  reconstructed  field  (trapezoidal  truncation,  M=25, 
L=25).  Let' 8 denote  its  measure  of  error  by 


AJ  l-^l  >2-5 


= Jk  \ j ( x - Xtf'wecM6 


an  unknown  quantity. 
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Magnitudes  of  Test  Field  Expansion  Coefficients  (xlO) 


Figures  9-4  illustrate  the  reconstruction  with  the  truncation  used  for 
the  predictions  made  to  date.  Trapezoidal,  M=9,  L=9 . The  error,  in 
terms  of  true  coefficient  values  including  negative  m's  is 
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6 - G-U'-is  + 50^4. 


Figures  9-5  thru  9-10  show  some  reconstructions  with  varying  truncation 


parameters . 

They  are  summarized 

below. 

Figures 

Truncation 
Policy  M 

L_ 

or  J 

Equivalent 
J or  L 

^ increment 

9-5 

Trapezoidal  9 

L 

= 15 

J 

= 10.5 

2122. 

9-6 

Rhomboidal  9 

J 

= 9 

L 

« 14. 

3175. 

9-7 

Trapezoidal  15 

L 

= 15 

J 

- 7.5 

1578. 

9-8 

Rhomboidal  15 

J 

= 9 

L 

* 17. 

2328. 

9-9 

Trapezoidal  15 

L 

= 25 

J 

- 17.5 

88. 

9-10 

Trapezoidal  12 

L 

o 

CM 

H 

J 

= 13.5 

623. 

i 


The  error  measures  for  these  reconstructed  fields  are  equal  to  ^-zc;  iS 
G increment  from  the  table. 

We  feel  that  the  case  illustrated  in  Figure  9-10  is  a reasonable  policy 
to  try  on  our  next  series  of  runs  (trapezoidal,  M=12,  L=20) . The 
G increment  for  the  equivalent  rhomb iodal  policy  is 


985.  , (Rhomboidal,  M-12,  >13) 

or  730.  , (Rhomboidal,  M*12,  >14) 

which  are  similar  to  the  G increment  (623.)  for  the  trapezoidal  case. 


Since  our  programs  are  not  set  up  for  rhomboidal  truncation,  we  have 
pursued  the  rhombiodal/trapezoidal  comparison  only  to  the  extent  of 
demonstrating  that  our  trapezoidal  policy  is  a reasonable  alternative. 
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(b) 


Figure  9-5 

(a)  Test  Field,  Trapezoidal  Truncation  ( M=  9,  L=15) 

(b)  With  Untruncated  Map  Overlay 
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Test  Field,  Rhomboidal  Truncation  ( 9 

With  Untruncated  Map  Overlay 


(b) 


Figure  9-8 

(a)  Test  Field,  Rhomboidal  Truncation  (M-15,  J=  9 ) 

(b)  With  Untruncated  Map  Overlay 
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Figure  9-10 

(a)  Test  Field , Trapezoidal  Truncation  ( M=  12,  I, = 20) 

(b)  With  Untruncated  Map  Overlay 
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Section  10.  TIKE  STEP  EXPERIMENT 


The  following  four  maps  represent  the  results  of  an  experiment  to  determine 
whether  an  increase  in  the  time  step  will  make  the  computation  blow 
up.  Four  twelve  hour  predictions  were  made.  The  time  step  used  were 
10  minutes,  20  minutes,  30  minutes  and  40  minutes.  We  display  the  results 
for  the  500  mb . pressure  height.  As  one  can  see  the  first  three  maps 
are  almost  identical.  The  40  minutes  time  step  differs  considerable 
from  the  other  three.  As  far  as  the  expansion  coefficients  are  con- 
cerned, the  coefficients  for  the  first  three  runs  differ  very  little  from 
each  other.  For  the  40  minute  runs,  the  low  M and  L coefficients  are  very 
similar  to  those  of  the  first  three  predictions,  those  with  higher  M and 
L values,  above  5 say,  are  considerably  different.  Values  of  all  co- 
efficients and  their  time-derivatives  have  been  checked  at  all  layers  as 
has  the  value  of  the  divergence  of  the  velocity  on  gridpoints  in  all 
layers.  There  Is  no  indications  of  any  instability  setting  in.  We 
thus  think  that  the  difference  as  far  as  the  40  minute  time  step  is 
concerned  is  simply  a matter  of  an  inaccuracy  in  integration  produced 
by  the  time  step  being  too  long. 
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500  mb,  00Z  2 Dec.  1974,  12  Hour  Prediction  (10  min.  time  step) 


500  mb,  00Z  2 Dec.  1974,  12  Hour  Prediction  (20  min.  time  step) 


500  mb.  00Z  2 Dec.  1974,  12  Hour  Prediction  ( 40  min.  time  step) 


Section  11.  ANGULAR  MOMENTUM 


The  angular  momentum  vector  of  a moving  fluid  of  density  p and  velocity 
V is  given  by  * 


where  d^is  the  volume  element  and  r is  the  radius  vector.  To  begin  with 
consider  the  case  of  a sphere  surrounded  by  a atmosphere  of  spherically 
symmetric  density.  '1  can  be  expanded  in  the  form 
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The  contributions  for  these  three  kind  of  terms  are  then 
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Consider  the  integral 
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Using  the  orthonormality  of  the  Y^'s  the  only  non-zero  terms  for  this 
integral  must  be  those  originating  in  C^’s  which  contain  y[J  terms.  A 
look  at  equation  (8)  shows  that  there  are  no  such  terms. 
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Consequently  LR  ■ 0. 
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Again  only  terms  In  containing  will  contribute.  These  terms  are 


99 


(13) 


~\ 


. \ 

' C 


r 


/\ 


r 

* V 

l c e 


V 


' \ \ c 


e 


or 


& 


i 


\ 


' c 


c 

\0 


\ 

' 0 


(14) 


Thus  we  have  reached  the  result  that  to  the  approximation  that  (■  i ‘ * i 
the  only  terms  capable  of  having  a total  angular  momentum  different  from 
zt»rc  are: 
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Now,  of  course,  two  Important  approximations  have  been  made:  the  earth 
was  assumed  spherical  and  the  density  was  assumed  to  depend  only  on 
the  distance  from  the  center  of  the  earth  but  not  on  the  angles  and  . 
These  are,  however,  pretty  good  approximations  and  to  a considerable  degree, 
then,  the  angular  momentum  of  the  atmosphere  resides  in  the  three  terms 
found  above.  In  particular  the  term  c^(r,t)  C^(fc,  ^)  is  the  main  component 
of  the  zonal  circulation. 
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The  imaginary  value  does  not  matter,  it  simply  means  that  c^(r,t)  also 
has  to  be  imaginary  to  make  the  velocity  and  the  angular  momentum  real. 

The  result  implies  that  must  be  very  stable  since  if  we  assume  that 
at  one  instant  of  time  the  angular  momentum  only  has  a z-component, 
i.e.  c^  (r,t)  = c ^ (r,t)  ■ 0,  and  since  the  other  velocity  terms  cannot 
absorb  any  total  angular  momentum,  the  law  of  conservation  angular  momentum 
will  keep  going.  We  show  below  some  results  from  real  data  and 
predictions 
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Section  12.  SELECTED  EXPANSION  COEFFICIENTS 


Out  of  the  very  large  number  of  expansion  coefficients  that  are  contained 
in  the  computer  printouts  we  have  chosed  to  display  graphically  a few 
generated  from  25  July,  1974  NMC  data.  This  is  so  that  the  reader  can 
get  a feel  for  how  they  depend  on  M and  L.  For  the  graphs  only  M is 
displayed,  the  value  of  L in  each  group  with  a given  M value  goes  from 
M to  9.  Of  the  wind  components  only  the  C®  coefficients  (graph  c“ 
expansion  coefficient  magnitudes)  are  displayed.  C°  represent  the 


zonal  components  and  C^,  (Mj*0)  represent  Rossby,  Haurwitz,  or  Neamtam 


waves,  whatever  you  choose  to  be  the  proper  nomenclature. 
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Section  13.  RESULTS  FROM 


10  MINUTES  PREDICTIONS 


We  present  in  the  following  pages  a few  selected  maps  of  the  500  mb 
temperature.  As  has  been  discussed  before  these  results  are  with  trun- 
cation limits  L = M ■ 9.  As  the  previous  discussion  on  truncation 
max.  max. 

limits  shows  that  these  are  not  reasonable  truncation  limits.  However, 
in  order  to  save  computer  time  we  decided  that  it  was  best  to  first 
ascertain  with  these  truncation  limits  that  the  time  step  was  reasonable, 
that  the  model  was  stable  and  that  the  results  looked  reasonable.  The 
model  has  been  run  out  to  5 days  or  144  hours  or  844  time  steps.  There 
seems  to  be  quite  a good  agreement  between  truncated  data  and  prediction 
for  pressure  height  to  48  hours  and  for  temperature  for  a few  days  more 
The  CPU  time  on  1MB  360/195  is  about  18  minutes  for  a 24  hour  prediction 
with  10  minutes  time  step  with  code  that  has  not  been  optimized  for  speed 
but  for  ease  of  change  and  flexibility.  We  are  using  the  transform  method 
from  coefficient  to  physical  space  at  every  time  step  but  not  the  Fast 
Fourier  Transform  since  we  did  not  want  to  be  tied  down  in  our 
experiments  to  the  gridspacing  required.  Provisions  have  been  made  in 
the  program  for  the  inclusion  of  the  F.F.T.  No  particular  initialization 
has  been  made  for  our  purposes.  We  have  simply  used  the  N.M.C  tapes 
in  the  form  they  have  been  used  by  N.M.C.  No  oscillations  to  speak  of  seem 
to  occur  in  the  first  12  hours,  which  is  where  we  have  monitored  selected 
data  for  every  time  step.  We  have  computed  the  kinetic  energy  of  the  zonal 
components  every  six  hours  and  the  fluctuation  is  of  the  order  of  0.5%. 

We  have  monitored  the  spherically  symmetric  part  of  the  logarithm  of  the 
surface  pressure  which  of  course  is  a measure  of  the  total  weight  of  the 
atmosphere  and  found  the  r.m.s.  deviation  from  the  mean  to  be  less  than 
one  part  in  a million  for  5 days.  We  have  at  every  six  hours  printed  1 out 
the  value  of  the  fivergence  of  the  velocity  at  every  gridpoint  in  every 
layer  and  found  no  significant  oscillations  except  that  the  predicted 
values  differ  from  the  data  about  one  order  of  magnitude.  The  values  in 
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the  data  always  being  smaller.  Apparently  the  divergence  is  being 
suppressed  in  the  data.  We  just  after  these  few  comments  and  maps  at 
this  time  in  order  to  keep  this  report  to  reasonable  lengths.  We  intend 
to  publish  a more  detailed  report  after  we  have  had  more  time  to  analyze 
both  the  data  and  the  prediction  and  after  we  have  extended  the  truncation 
limits. 

The  maps  represent  the  500  mb  pressure  height.  NMCIC  represents' a map 

printed  out  from  N.M.C  is  data  tape.  Prediction  is  the  same  pressure 

height  with  L = M =9  and  NMC  I.C.  truncated  means  that  the  data 
max.  max. 

has  been  converted  to  expansion  coefficients  for  the  spherical  harmonics 
with  the  same  truncated  limits  and  then  converted  to  numerical  values 
from  the  functions  on  the  grid. 
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